Exact measurement of the properties of nature is a common goal within the experimental sciences. Similarly, medical diagnostics and communications technology, among other scientific endeavors, seek the ability to obtain exact measurement of properties within their respective fields, e.g., MRI or free space optical or other electromagnetic transmission. Optimal extraction of the underlying measurement data requires the removal of measurement defects such as noise and limited instrumental resolution. However, in spite of the availability of highly sophisticated instruments, noise and instrumental signatures are still present in the data, making the measurement only approximate.
An area of experimental science in which instrumental signatures are particularly a problem is astronomy, where the signal sources to be measured are very faint. Even when the instruments are made essentially noise-free, instrumental signatures related to finite spatial, spectral, or temporal resolution remain. At this point, image reconstruction is required to remove the instrumental signatures.
One of the most powerful approaches to image restoration, e.g., removal of point-spread-function blurring, is Bayesian image reconstruction, which includes Goodness-of-Fit (Maximum Likelihood) and Maximum Entropy data fitting. This family of techniques employs a statistical relationship between various quantities involved in the imaging process. Specifically, the data, D, consisting of the original noisy, blurred image model is linked to the desired noise-free, unblurred image model, I, through a model, M. The model M includes all aspects of the relationship between the data and the image, e.g., that the data is normally collected on a rectangular grid and that the data is related to the image model through the relationshipD(i)=(I*H)(i)+N(i)  (1)where D(i) is the data in cell i (typically a pixel), I is the image model, H is the point-spread-function (PSF) due to instrumental and possible atmospheric blurring, * is the spatial convolution operator, i.e.,                                           (                          f              *              g                        )                    ⁢                      (            x            )                          =                              ∫                          -              ∞                        ∞                    ⁢                                           ⁢                                    ⅆ                              x                ′                                      ⁢                          f              ⁡                              (                                  x                  ′                                )                                      ⁢                          g              ⁡                              (                                  x                  -                                      x                    ′                                                  )                                                                        (        2        )            and N represents the noise in the data, assuming the PSF is a function only of displacement between pixels. In general, the PSF can vary across the field.
Image reconstruction differs from standard solutions of integral equations due to the noise term, N, that nature of which is only known statistically. Methods for solving such an equation fall under the categories of (1) direct methods, which apply explicit operators to data to provide estimates of the image, and (2) indirect methods, which model the noiseless image, transform it only forward to provide a noise-free data model, then fit the parameters of the image to minimize the residuals between the real data and the noise-free data model. The direct methods have the advantage of speed, but they tend to amplify noise, particularly at high frequencies. The indirect methods supposedly exclude the noise, however, the required modeling is a disadvantage. If a good parametric form for the image is known a priori, the result can be very good.
To statistically model the imaging process, Bayesian image reconstruction methods analyze the properties of the joint probability distribution of the triplet, D, I and M, i.e., p(D,I,M). Applying Bayes' Theorem [p(A,B)=p(A|B)p(B)=p(B|A)p(A), where p(X|Y) is the probability of X given that Y is known] provides:                                                                                                                                                                                       p                          ⁡                                                      (                                                          D                              ,                              I                              ,                              M                                                        )                                                                          =                                                ⁢                                                                              p                            (                            D                                                                                ⁢                          I                                                                    ,                      M                                        )                                    ⁢                                      p                    (                    I                                                        ⁢                                                                           ⁢                  M                                )                            ⁢                              p                ⁡                                  (                  M                  )                                                                                                                                                                                                            =                                                ⁢                                                                              p                            (                            I                                                                                ⁢                          D                                                                    ,                      M                                        )                                    ⁢                                      p                    (                    D                                                        ⁢                  M                                )                            ⁢                              p                ⁡                                  (                  M                  )                                                                                        (        3        )            By setting the first factorization of p(D,I,M) in Equation 3 equal to the second factorization provides the usual starting point for Bayesian reconstruction:                                                                         p                (                I                                            ⁢              D                        ,            M                    )                =                                                                                                                        p                      (                      D                                                              ⁢                    I                                    ,                  M                                )                            ⁢                              p                (                I                                            ⁢              M                        )                                                              p                (                D                                            ⁢              M                        )                                              (        4        )            
A common goal of Bayesian image reconstruction is to find the M.A.P. (Maximum A Posteriori) image, I, which maximizes p(I|D,M), i.e., the most probable image given the data and model. (Note that other image estimates, e.g., the average image, <I>=∫D,M dMdD I p(I|D,M), may be used here and in the methods described in the detailed description.) (MAP image reconstruction is also sometimes known as parametric least-squares fitting.)
It is common in Bayesian image reconstruction to assume that the model is fixed. In this case, p(D|M) is constant, so thatp(I|D,M)∝p(D|I,M)p(I|M).  (5)
The first term, p(D|I,M), is a goodness-of-fit quantity, measuring the likelihood of the data given a particular image and model. The second term, p(I|M), is normally referred to as the “image prior”, and expresses the a priori probability of a particular realization of the image given the model. In Goodness-of-Fit (GOF) image reconstruction, p(I|M) is effectively set to unity, i.e., there is no prior bias concerning the image. Only the Goodness-of-Fit (p(D|I,M)) is maximized during image reconstruction. Typically,p(I|D,M)=exp(−χ2R/2),  (6)where χ2R is the chi-square of the residuals, R(≡D−I*H). While this approach ensures that the frequency distribution of the residuals has a width which is characteristic of the noise distribution, it normally results in images with spurious spatial features where the data has a low signal to noise ratio (SNR). Also, the large amplitude residuals often show a strong spatial correlation with bright features in the data.
When no parametric model of the image is known, the number of image model parameters can quickly become comparable to, or exceed, the number of data points. In this case, a MAP solution becomes inappropriate. For example, if the number of points in the image model equals the number of data points, the nonsingular nature of the transform in Equation 1 assures that there is a solution for which the data, including the noise, are exactly modeled with zero residuals. This is the same poor solution, with all its noise amplification, obtained by the naive Fourier deconvolution. Thus, an unrestricted indirect method is no better at controlling noise than a direct method, and therefore, the image model must be restricted in some way. The indirect methods restrict the image model and differ only in the specifics of image restriction.
A simple restriction is to constrain the model image to be positive. Since even a delta-function image is broadened by the PSF, it follows that the exact inverse of any noisy data with fluctuations on scales smaller than the PSF must be both positive and negative. By preventing the image from becoming negative, the noise-free data model cannot fluctuate on scales smaller than the PSF, which is equivalent to smoothing the data on the scale of the PSF. However, this Non-Negative Least-Squares (NNLS) fit method is not able to eliminate noise fitting on larger scales.
Maximum entropy (ME) image reconstruction solves many of the problems of the simpler GOF methods (e.g., NNLS). In ME imaging, one calculates a value for the image prior based upon “phase space volume” or counting arguments. Heuristically, p(I M) is written p(I M)=exp(S), where S is the entropy of the image in a given model. All ME methods capitalize on the virtues of incorporating prior knowledge of the likelihood of the image. The information entropy is maximized for a flat image, which eliminates structure not required by the data and suppresses noise fitting. The benefits of this are numerous, including eliminating the over-resolution problems of GOF methods and increasing the numerical stability of the calculations. The disadvantages are that the image prior used in ME takes the form of a global constraint, resulting in the enforcement of an average smoothness (average information content) on the entire image without recognizing that the density of information content, and, thus, the appropriate degree of smoothing, varies from location to location in the image. This results in oversmoothing of the image in some locations and undersmoothing it in others.
Many Bayesian image reconstruction methods assume that the model is fixed. Others, however, propose varying the model, such as the multi-channel image modeling of Weir (Applications of Maximum Entropy Techniques to HST Data, Proceedings of the ESO/ST-ECF Data Analysis Workshop, April 1991). In this method, the image is assumed to be a sum of pseudoimages convolved with a blurring function of various spatial scales. This method, while superior to many of its predecessors, may exhibit low-level spurious sources as ripples in the reconstructed image, and still displays some spatial correlation within the residuals.
The need remains for a method of image reconstruction which is capable of effectively reducing noise fitting without relying on a priori probabilities that result in excessive averaging of the image, but rather which adapts itself to the distribution of information content in the image. The following disclosure provides such a method.